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Abstract. We consider the slow nonlinear diffusion equation subject to a constant absorption rate and construct local 
self-similar solutions for reversing (and anti-reversing) interfaces, where an initially advancing (receding) interface gives way 
to a receding (advancing) one. We use an approach based on invariant manifolds, which allows us to determine the required 
asymptotic behaviour for small and large values of the concentration. We then ‘connect’ the requisite asymptotic behaviours 
using a robust and accurate numerical scheme. By doing so, we are able to furnish a rich set of self-similar solutions for both 
reversing and anti-reversing interfaces. 
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1. Introduction. We address reversing and anti-reversing properties of interfaces in the following 
one-dimensional slow diffusion equation with strong absorption 


^ (h^—\ - 

dt dx I dx' 


( 1 . 1 ) 


where h is a positive function, e.g., a concentration of some species, and x and t denote space and time, 
respectively. Restricting the exponents to the ranges m > 0 and n < 1 corresponds to the slow diffusion 
and strong absorption cases respectively. 

Interfaces — sometimes termed ‘contact lines’ by fluid dynamicists — correspond to the points on 
the x-axis, where regions for positive solutions for h are connected with the regions where h is identically 
zero. The initial data h\t=o = Hq is assumed to be compactly supported. The motion of the interfaces 
is determined from conditions that require the function h be continuous and the flux of h through the 
interface to be zero [5]. 

In the presence of slow diffusion (m > 0), the interfaces of compactly supported solutions have a finite 
propagation speed m- In the presence of strong absorption (n < 1), the solution vanishes for all x after 
some finite time, which is referred to as finite-time extinction [MlIIl]. Therefore, the interfaces for a 
compactly supported initial data coalesce in a hnite time. Depending on the shape of ho and the values of 
m and n, the interfaces may change their direction of propagation in a number of different ways. It was 
proved by Chen et al. [3] that bell-shaped initial data remains bell-shaped for all time before the compact 
support shrinks to a point. However, the possible types of dynamics of interfaces for this bell-shaped data 
were not identified in [3]. 


The slow diffusion equation with the strong absorption (1.1) describes a variety of different physical 


processes, including: (i) the slow spreading of a slender viscous film over a horizontal plate subject to the 
action of gravity and a constant evaporation rate ^ (when m = 3 and n = 0); (ii) the dispersion of a 
biological population subject to a constant death-rate |10j (when m = 2 and n = 0); (iii) non-linear heat 
conduction along a rod with a constant rate of heat loss [11] (when to = 4 and n = 0), and; (iv) fluid 
flows in porous media with a drainage rate driven by gravity or background flows [5J [T^ (when to = 1 and 
either n = 1 or n = 0). 

Let us denote the location of the left interface by x = £{t) and the limit x \ £{t), where h is nonzero, 
by X = £{t)^- If TO -I- n > 1, it was proved in [^ that the position of the interface, £{t), is a Lipschitz 
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continuous function of time t. In the case m + n = 1, the function £(t) is found from the boundary 
conditions = 0 and 


OX 




,m-l^ 

dx 




( 1 . 2 ) 


where a dot denotes differentiation with respect to time. In the case to + n > 1, the spatial derivatives at 


X = £{t)'^ are not well defined [5]. and the zero flux condition (1.2) must be rewritten as 


i = 


-h' 


h'^ 


^-idh 

dx 


dh 

dx 


c=l{ty 




if £ < 0, 


if £ > 0. 


(1.3) 


One could choose to close the slow diffusion equation (O) in a variety of ways, e.g., by supplying 
analogous conditions at the right interface, or by supplying a Dirichlet or Neumann condition elsewhere. 
For instance, if /iq is even in x, then the solution h remains even in x for all times, and therefore, the 


slow diffusion equation (1.1) can be closed on the compact interval [£(£), 0] by imposing dhldx\x=o = 0. 
However, such details do not concern us here because we are interested in studying the behaviour of 


solutions to (1.1) local to the left interface x = i{t) only. 


We reiterate here the main question on the possible types of dynamics in the slow diffusion equation 


with the strong absorption (1.1). Working with bell-shaped, compactly supported initial data Hq, one 


can anticipate a priori that the compact support of the bell-shaped solution can either: (i) decrease 
monotonically in time, or; (ii) first expand and then subsequently shrink, or; (iii) have more complicated 
behaviour where multiple instances of expansion and contraction are observed. This phenomenon brings 
about both ‘reversing’ and ‘anti-reversing’ dynamics of an interface. Here the term ‘reversing’ describes a 
scenario where the velocity of the left interface x = £{t) satisfies £ < 0 before the reversing time and £ > 0 
after the reversing time, whereas the term ‘anti-reversing’ refers to the opposite scenario with £ > 0 before 
and £ < 0 after the reversing time. 


The first analytical solution to (1.1) exhibiting a reversing interface was obtained by Kersner [T3] for 
the case m + n = 1. This explicit solution takes the form 


\x,t) = 


2(to -|- 1)(to -|- 2)t - 


Ct^+^ - {m + 2)^t^ - x"^ 


J -t 


(1.4) 


where the plus subscript denotes the positive part of the function, and C > 0 is an arbitrary parameter. 
The interfaces are located symmetrically at a: = ±£(£) with 


£(£) = Ct^+^ — (m + 2yt'^. 


(1.5) 


More recently, Foster et al. [B] considered the case to -I- n > 1 and explored the asymptotic and 
numerical construction of self-similar solutions for equation (1.11 — some related, yet different, self-similar 


solutions to other nonlinear diffusion equations have previously been constructed using a combination of 
asymptotic analysis and numerical shooting; see, e.g., laiiH!- 

The self-similar solutions capture the relevant dynamics of reversing interfaces near the corresponding 
points in the space and time (which can be placed at the origin of x and £, without the loss of generality). 


Based on a classical point symmetry analysis of the porous medium equation (1.1) — provided in — 
the authors of [B] found that the reversing interfaces can be described via the self-similar reductions 


h{x,t) = ^ = x{±t) 2 ( 1 -™)^ ±i > 0, 


( 1 . 6 ) 
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where the functions H± satisfy a pair of second-order differential equations. 

In this paper, we explore the case n = 0 only, when the functions H± of the self-similar reduction 


(1.6) satisfy the second-order differential equations 




dH±\ m + 1 dH± 


d^ 


± 


d^ 


= 1±H±. 


(1.7) 


We seek positive solutions H± of the differential equations (1.7) on the semi-infinite line [A±,oo) that 
satisfy the following conditions: 


(i) 

(ii) 

(iii) 


77±(^)—as ^ t 7l±, 

H±{0 is monotonically increasing for all ^ > A±, 
^ -foo as ^ —y -t“Oo. 


( 1 . 8 ) 

(1.9) 

( 1 . 10 ) 


These first of these conditions, (1.8), are the self-similar counterparts of the condition = 0 for the 

equation (1.1). In addition, the behaviour of 77-(C) and in the far-held must be matched from the 

condition 


H+iO 

C->oo H_{^) 


lim 


= 1 . 


( 1 . 11 ) 


The requirement (1.11) is tantamount to enforcing that the solution h to the slow diffusion equation (1.1) 


does not ‘jump’ as t passes through zero — this can be verihed by taking both the limits t \ 0 and t Z' 0 


in the self-similar reduction (1.6). 


Existence of solutions to the differential equations (1.7) on [A±, oo) with the required behaviour (1.11) 
implies, via the self-similar reduction (1.6), the existence of a reversing (if A± > 0) or anti-reversing (if 
A± < 0) left interface at cc = i{t), which behaves like 


£{t) = A±{±t)^ 


> 0 , 


( 1 . 12 ) 


after placing the reversing point at the origin of the space-time plane. If m > 1, the velocity of the 
interface, £{t), changes sign continuously as t passes through zero. 

By combining formal asymptotic constructions of the solutions H± near the small and large values 
with a numerical shooting method, the authors of [S] claimed that for integer values of to = 2,3,4, there 
exists a unique positive value of A-, which leads to a monotonically growing function i7_ on the entire 


semi-axis [A-, oo). Using the far-held matching condition (1.11), a unique, monotonically growing function 
i7_i_ is found on [A_|_, oo) for a positive value of A+. These solutions correspond to a reversing left interface 
X = £(t) for TO = 2,3,4. No solutions exhibiting an anti-reversing interface were found in [5]. 

In the present work, we address the same problem using a dynamical system framework dUn]. The 
dynamical system theory allows us to justify the formal asymptotic approximations of H± for small and 
large values, as well as to set up an accurate and robust numerical procedure for furnishing appropriate 


solutions to the differential equations (1.7). Qualitatively, we recover the results of [5] for to = 3,4, but 
with a better precision, and we generalize these results for all non-integer values of to > 1. In addition, we 
demonstrate that the result for to = 2, reported in [^, is incorrect and no self-similar reversing interface 
solutions exist for to = 2. In addition, we discover new reversing and anti-reversing interface solutions of 


the same differential equations (1.7) for other values of to. 


Our approach explores invariant manifolds for the singular differential equations after appropriate 
unfolding (which is sometimes referred to as the blow-up technique [51 116jl. The main analytical results 
of this work are given by the following theorems. 


Theorem 1.1. For every to > 1 and A± ^ 0, there exists a unique solution of the differential equation 
0 such that H±{f) —>■ 0 as ^ A±. If ±A± > 0, this unique solution has the following asymptotic 
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behaviour 




(C-A±) + 0 ((e-A±r“{ 2 .m})^ as 


' (m + l)A± 

whereas if iA± < 0, it has the following asymptotic behaviour 


(1.13) 


(1.14) 


Theorem 1.2. There exists a one-parameter family of solutions of the differential equation (1.1) for 
the lower sign such that i7_(^) —> +oo as ^ > +oo, and this family has the following asymptotic behaviour 


H-iO = - 

\xo 


(l + C>(^ ^)) , as ^+ 00 , 


(1.15) 


where Xq > 0 is an arbitrary parameter. There exists a two-parameter family of solutions of the differential 
equation (1.1) for the upper sign such that T{+{f) —)■ +c» as ^ +oo, and this family has the same 


asymptotic behaviour (1.15) for some Xq > 0. 


The main problem is to connect the two asymptotic behaviours of the differential equations (1.7) 
which are defined for small and large values of H± by Theorems and |1.2[ We know from [6] that there 
exists an exact solution of the connection problem if A± = 0. This exact solution is given by 


H ±{0 = 


m + 1 




f. e (0,oo). 


(1.16) 


However, the exact profile (1.16) corresponds to a solution to the slow diffusion equation (1.1) with an 


interface that remains stationary for all time, thus we do not examine it further here. 

Some connection results for nonzero values of A± are available for the differential equation (1.7) in 
Lemmas 14.11 and 14.41 below. 


Although these results are not sufficient for an analytical solution of the 
connection problem, we can set up a numerical method, which detects connections of the two solutions 
described in Theorems 11.11 and 11.21 

The remainder of the paper is organized as follows. The unfolding and invariant manifolds for the 


differential equations (1.7) near small values of H± are described in §2. The corresponding results near 
large values of H± are reported in §3. The connection problem between the invariant manifolds near small 
and large values of H± is considered in §4. The relevant numerical technique is implemented in §5, where 
the main findings are discussed and compared with the previous results from [5]. The paper is concluded 
in §6 with a discussion of the relevance of the self-similar solutions to the dynamics of (1.1). 


2. Invariant manifolds for small values of H±. We shall rewrite the scalar equations (1.7) as 
vector systems for variables u = H± and w = Hlff . In the interests of simplicity of notation, we drop 
the plus and minus subscripts in the definitions of the variables u and w. The non-autonomous vector 
system for u and w is as follows: 


du _ w 

= 1±UT 

d^ ' 2 


( 2 . 1 ) 


If TO is a non-integer, we require the constraint u > 0. In either case, only positive solutions for u are 
needed to be considered. 

The interface, in self-similar variables, is assumed to be located at ^ = A G K, where u = 0 — a 
requirement of the condition on the continuity of h{x, t) ai x = £(t). Since the value u = 0 is singular in the 
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non-autonomous system ( 2 . 11 , we shall unfold the singularity by introducing a convenient parametrization 


of solutions with a new time variable r defined by 


dr 


= u 


u> 0. 


The map r i—^ is increasing and if ^ > A as r —> —oo, then ^ ^ for finite values of r. 

With the parametrization r i-A we obtain the autonomous dynamical system in 


ii = w, 


( 2 . 2 ) 


w 


where the dots stand for the derivatives of {^,u,w) in r. In what follows, we assume that m > 1, so that 


the vector field of the dynamical system (2.2) is continuously differentiable near zero values of u. Again, 


u > 0 has to be enforced if m is a non-integer. 

The family of equilibrium points for the system (2.2) is given by (^,u,w) = (A, 0,0), where A € 
an arbitrary parameter. If m > 1, each equilibrium point is associated with the Jacobian matrix 


IS 


0 0 0 
0 0 1 
0 0 t^A 


This Jacobian matrix has a double zero eigenvalue (with two linearly independent eigenvec tors ) and a 
simple nonzero eigenvalue Therefore, the linearization of the dynamical system (2.2) at the 

equilibrium point (A, 0, 0) with ±A > 0 has a two-dimensional center manifold and a one-dimensional 
stable manifold, whereas the linearized system with ±A < 0 has a two-dimensional center manifold and 
a one-dimensional unstable manifold. Since the dynamical system (2.2) is smooth, a straightforward 
application of the invariant manifold theorems da Hi] asserts that the equilibrium point (A, 0,0) with 
A 7 ^: 0 is located at the intersection of the two invariant manifolds, which are tangential to the invariant 
manifolds of the linearized system. We formulate these results in the two following Propositions, namely 
Propositions |2.1| and |2.2| The relevant conclusions on the behaviour of solutions of the differential equa¬ 
tions (1.7) for small values of ff±, expressed in Theorem [Hi follow from these two Propositions. 


Proposition 2.1. For every m > 1 and A 7^ 0, there exists a two-dimensional center manifold of the 
dynamical system (2.2) near the equilibrium point (A, 0,0), which can be parameterized as follows: 


W,(A,0,0) = \w = ± 


2m” 


' (m -I- 1)A . 


1 + 0(C-A,m““{1’'"-i>) 1 , MS (0,(5), ^e(A,A + (5), 


(2.3) 


where 5 > 0 is small. Dynamics of the system (2.2) on the center manifold Wc(A,0,0) is topologically 
equivalent to the dynamics at the truncated normal form 


f = u^, 
u = ± 


(m+l)A ' 


(2.4) 


In particular, for every A 7 ^ 0, there exists exactly one trajectory on Wc(A, 0,0), which approaches the 
equilibrium point (A, 0,0) as t ^ —00 if ±A > 0 and t —>■ J-cc if ±A < 0. 

Proof. Existence of a two-dimensional center manifold ITc(A, 0, 0), which is tangent to that of the 
linearized system 


Ec(A, 0,0) = {w = 0, (^,m)gM^}, 
















6 


J. M. FOSTER AND D. E. PELINOVSKY 


follows from Theorem 4.1 in [3]. We develop an approximation of Wc(A,0,0) by writing 

w = 


(2.5) 


where u) is at the point (^, u) = {A, 0) with zero partial derivatives. Dynamics along Wc{A, 0, 0) 

is given by the two-dimensional system 


ii = u). 


( 2 . 6 ) 


two-dimensional system (2.6). Then, we obtain a partial differential equation 


The function ry is to be found by substituting (2.5) to the three-dimensional system (2.2) and using the 

(2.7) 


It ^^Ar] = Tu±^^{^- A)r] + r]—{u ??)-h — (u ry). 


If m > 1 and u) is a function at (^, u) = {A, 0) with zero partial derivatives, then equation (2.7) 

has a solution such that 


u)=± 


{m -I- 1)A 


-FC>(C- 


( 2 . 8 ) 


The representation (2.5) and (2.8) is equivalent to (2.3). Substituting (2.8) to (2.6) and truncating the 
remainder term, we obtain the truncated normal form (2.4). 

From the second equation of the system (2.4), it follows that if zLA > 0, then u > 0 such that m(t) —>■ 0 
as r —>■ — 00 , wh ereas if ±A < 0, then u < 0 such that u{t) —)■ 0 as r —> -l-oo. From the first equation of 
the system (2.4), the constant of integration for ^ is arbitrary, so that ^(t) A in the same limit with 
A^ A. Hence, dynamics along the two-dimensional manifold Wc(H, 0,0) is decomposed between a curve 
of equilibrium states with u = 0 and weakly unstable (if zLA > 0) or weakly stable (if zLA < 0) evolution 
along a curve parameterized by small positive u. 

Persistence of the dynamics on Wc(H, 0, 0) with respect to the remainder terms in (2.8) follows from 
analysis of the system (2.6). □ 


Proposition 2.2. For every m > I and zLA < 0, there exists a one-dimensional unstable manifold 

which can be paramet 

+ 0(u™), uG{0,S), 


of the dynamical system [2.2) near the equilibrium point (H,0, 0), which can be parameterized as follows: 

2u 


Wu{A,Q,Q) = \i = A + 0{u^), w = T 


(to -|- 1)A 


(2.9) 


where <5 > 0 is small. Dynamics of the system 12.^ on the unstable manifold W„(H,0,0) is topologically 
equivalent to dynamics of the linear equation 


TO -I- 1 

u = —-— Au. 


( 2 . 10 ) 


Proof. Existence of a one-dimensional unstable manifold IF„(H, 0,0), which is tangent to that of the 
linearized system 

I jYi -j- 1 

Etj(H,0,0) = <^ ^ = H, w = T —X— Au, uG 


follows from Theorem 4.1 in [4]. We develop an approximation of Wu{A,0,0) by writing 

^ H -|- 


_ -P m+1 


W = ^ 


Au u^9{u), 


( 2 . 11 ) 
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where u'^4>{u) and u'^6{u) are with the zero derivative at m = 0. Dynamics along 0, 0) is given 

by the one-dimensional system 


W + 1 < mn, \ 

u = =F—s —Au + u 9{u). 


( 2 . 12 ) 


using the one-dimensional system (2.12). Then, we obtain a system of differential equations 

d(j)^ 


The functions (j) and 9 are to be found by substituting (2.11) to the three-dimensional system (2.2) and 

ysten 

) = 1 


m<p{u) + u 


du 


m+1 ^ 


u 


(2.13) 


and 


T- 


1 


A- 


^9(u) 


n, ^ d9\ 

m9{u) + u— 
du J 


= 1±mH- —(t){u) 


m + 1 


Au =F u'^9{u) 


(2.14) 


If TO > 1 while u^(j){u) and u'^9{u) are with the zero derivative at rt = 0, then system (2.13) and (2.14) 
has a solution such that 


</>(«) = T- 


mim ■ 


i)A 


+ Oiu^-^), 9 {u) = t 


to(to -I- 1)A 


+ 0(u 


min{ l,m— 1} \ 


(2.15) 


The representation (2.11) and ( |2.15| ) is equivalent to (2.9). Substituting (2.15) for 9{u) to (2.12) and 
truncating the remainder terms at OivA^), we obtain the linear equation (2.10). Persistence of the linear 
dynamics on Wu{A, 0, 0) with respect to the remainder terms in 9{u) follows from analysis of the differential 
equation (2.12). □ 


Remark 1. For every to > 1 and FA > 0, one can construct a one-dimensional stable manifold 


II/s(^, 0,0) of the dynamical system (2.2) near the equilibrium point (A, 0,0), which exists for f > A, 
u > 0, and w < 0. However, this stable manifold does not contain trajectories that approach the equilibrium 
point (A, 0,0) as T ^ —oo. 


Proof of Theorem For every to > 1 and FA ^ 0, Propositio n |2.1| states that the equilibrium state 
(A, 0,0) is connected by the trajectories of the dynamical system (2.2) with u{t) > 0 as r —>■ —oo if and 
only if FA > 0. In this case, there exists exactly one trajectory with u > 0 such that u(r) —>■ 0 as r ——oo. 


This trajectory belongs to the center manifold IFc(^)0,0), whose dynamics satisfy the system (2.6). From 
this system, we obtain a first-order non-autonomous equation 


^ = Vitu) = ± 7 — 

df (to - 


1)^ 


FO{f- as 


0 . 


(2.16) 


Integrating (2.16) near ^ = A, we recover the asymptotic behaviour (1.13). 


For every to > 1 and FA < 0, Proposition 2.2 states that the equilibrium state {A, 0, 0) is connected by 
exactly one trajectory of the dynamical system (2.2) with u > 0 and u{t) —0 as t —>■ —oo. This trajectory 
belongs to the unstable manifold IF„(^, 0, 0) with the dynamics satisfying equation ( |2.12[ ). From (2.11) 
and (2.15), we obtain 


e = A + u” 


T 


to(to + \)A 


FO(u'^-^) 


as M —>■ 0. 


Inverting this nonlinear equation near f = A, we recover the asymptotic behaviour (1.14). 


(2.17) 

□ 



















































J. M. FOSTER AND D. E. PELINOVSKY 


3. Invariant manifolds for large values of H±. The trajectories departing from the equilibrium 
point (A, 0,0) of the dynamical system (2.2) is expected to arrive at infinite values for ^ and u. In order 
to study the behaviour of trajectories near infinite values for ^ and m, we shall define y = 1/u, which maps 
an infinite value for m to a zero value for y. The other variables ^ and w must be adjusted accordingly for 
small values of y. Let us consider the following scaling transformation, 


yP^ 


1 

U= 

y 


w = —, 


(3.1) 


where (x, y, w) is the set of new variables, and the positive parameters p and q are to be chosen below. 
Substituting the transformation (3.1) to the dynamical system (2.2), we obtain the following autonomous 
system in 


X = yP-"^ -pxzy^-'i, 
y = —zy‘^~^, 

i = 


■y)T '^xzy p - qz'^y^ «, 


(3.2) 


where a dot still denotes a derivative with respect to the time variable r. 

The system (3.2) is singular at ?/ = 0, no matter what positive values oip and q are chosen. To unfold 
the singularity, we can now introduce a convenient parametrization of solutions with the time variable s 
instead of the variable r such that 


dr 

ds 


= y^, 


y>0- 


The map s i-A r is increasing and we can consider solutions parameterized by the new time variable s such 
that 2 / —>■ 0 as s —>■ +oo (which could correspond to a finite value for the old time variable r). 

With the parametrization s i-A r, the system (3.2) can be rewritten in the equivalent form 


x'= y'^P ™ — pxzyP^^ 
y' = -zyP+'^-'i, 

z' = yP+9-™-i(±i + y) ip P^xz - qz^yP+^-<^, 


(3.3) 


where a prime now denotes a derivative with respect to the new time variable s. A suitable choice of 
parameters p and q is given by 


p+l = q, 
p + q = m + 2, 


_ m+l 
P — 2 ’ 

^ _ m+3 

H — 2 ’ 


(3.4) 


so that the explicit form of the transformations (3.1) is given by 




y~ 


1 

M = -, 

y 


w = 


(3.5) 


y~ 


The choice (3.4) ensures that the system (3.3) is rewritten in the simplest non-singular form with a 
quadratic vector field: 


P^XZ, 


X =y 

y' = -zy, 
z' = y{±l + y)T - p^z'^. 


(3.6) 


The family of equilibrium points for the system (3.6) is given by (x, y, z) = (xq, 0,0), where Xq S K is 
an arbitrary parameter. Each equilibrium point is associated with the Jacobian matrix 

0 1 


0 

±1 


2 ^0 
0 

T'^xo 
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This Jacobian matrix has a double zero eigenvalue (with two linearly independent eigenvectors) and a 
simple eigenvalue =F^^^a;o. 

In the context of reversing and anti-reversing interfaces, we are interested in the behaviour of solutions 
for which ^ —>■ -|-oo and u —>■ -foo. It follows from the transformation (3.5), that this requirement restricts 


our consideration to the family of critical points (a:o,0,0) with xq > 0. If xq > 0, the linearization of the 
dynamical system (3.6) at the equilibrium point (a;o,0,0) has a two-dimensional center manifold and a 


one-dimensional stable (upper sign) or unstable (lower sign) manifold. Since the vector field of the dy¬ 
namical system ( |3.6[ ) is analytic (quadratic), another straightforward application of the invariant manifold 
theorems [13 I17| yields that the equilibrium point (a;o,0,0) with a;o > 0 is located at the intersection 
of the two invariant manifolds, which are tangential to the invariant manifolds of the linearized system. 
We formulate these results in the next two Propositions, namely Propositions |3.1| and |3.2| The relevant 
conclusions on the behaviour of solutions of the differential equations (1.7) for large values of H± expressed 
in Theorem 1 1.2 1 follow as a corollary from these two Propositions. 


Proposition 3.1. For every xq > 0, there exists a two-dimensional center manifold of the dynamical 

I, 0), wi 

m -f 1 2 \ 2 


system {3.6} near the equilibrium point (a:o,0, 0), which can be parameterized as follows: 

m 1 


Wc(a;o,0,0) = <y = 


xz± { 1 — 


Xq] Z 0(3) 


X G {xo — S, Xo + J), 

(-J,J), 


(3.7) 


where J > 0 is small and 0(3) denotes cubic terms in x — xq and z. The dynamics of the system (3.6) on 
the center manifold Wc(aio,0, 0) is topologically equivalent to the dynamics at the truncated normal form 




z = —z. 


(3.8) 


In particular, there exists exactly one trajectory on Wc(a;o,0,0), which approaches the equilibrium point 
{xq, 0,0) as s ^ -foo. 

Proof. Existence of a two-dimensional center manifold Wc(xo,0, 0), which is tangent to that of the 
linearized system 


C — 1 — 

Ec{xo,0,0) = ly = —-— XqZ, (a:,2;)GM^ 

follows from Theorem 4.1 in We develop an approximation of lTc(a;o:0j0) by writing y = f{x,z) and 
expanding / using a Taylor series in small values of both x — xq and z. Dynamics along Wc{xq,Q,Q) is 
given by the two-dimensional system 


x' = f{x,z)-^xz, 

z' = fix, z)i±l + fix, z)) T ^XZ - ^z^. 


The function / is to be found from the partial differential equation 


dx 


m-\-l 

J - —xz 


dz 


TO-l-1 m-l-32 
/(±1 + /) T -^^xz -— 


+ zf = 0. 


(3.9) 


(3.10) 


Equations (3.9) and (3.10) suggest the following near-identity transformation of the function / given by 

fix,z) = '^^^ xz-\- z'^gix, z) (3.11) 






















10 


J. M. FOSTER AND D. E. PELINOVSKY 


After the near-identity transformation (3.11), the dynamics along Wc{xq,0^0) are given by the two- 
dimensional system 


c' = z'^9{x,z), 


Z = Z 


±g{x,z)-^ + {^+zg{x,z)y 


(3.12) 


The function g is now to be found from the partial differential equation 


m + 3 dq 
zg I —^-h z — 


1 


2 ' dxj ' \ 2 

which has a solution such that 


-X -I- 2zg + z — 


m -I- 1 



/ \ ' 
g[x,z) = ±- 


1 I 


x + zg] - 


xl + 0{x - xo,z). 


m + 3 


1 


-X = 0, 


(3.13) 


Expansions (3.11) and (3.13) yield (3.7). Substituting (3.13) into (3.12) and truncating at the cubic terms, 


we obtain the truncated normal form (3.8). 


From the second equation of the system (3.8), it follows that there exists a unique solution such that 


z(s) —0 as s —>■ - 1-00 with z > 0. From the first equation of the system (3.8), the constant of integration 
for X is arbitrary, so that a:(s) —>■ io as s —>■ -l-oo with xq ^ xq- Hence, dynamics along the two-dimensional 
manifold Wc(xo,0,0) is decomposed between a curve of equilibrium states with z = 0 and weakly stable 
evolution along a curve parameterized by small positive z. Persistence of the dynamics on Wc{xq,Q,Q) 
with respect to the remainder terms in (3.13) follows from analysis of the system (3.12). □ 


Proposition 3.2. For every xq > 0, there exists a one-dimensional stable (upper sign) or unstable 


(lower sign) manifold of the dynamical system (3.6) near the equilibrium point (xo,0,0), which can be 
expressed explicitly: 


I jYi 1 

W^s/«(a;o,0,0) = < y = 0, z = T—z —x 


1 - 


X € {xq — (5, xo + 6), 


(3.14) 


where <5 > 0 is small. The dynamics of the system (3.6) on the manifold HA/u(a;o,0,0) is topologically 
equivalent to the dynamics of the linear equation 


x' = T- 


m -|- 1 


xo{x - Xo). 


(3.15) 


Proof. The linearized system has the stable/unstable manifold: 

Es/uixo,0,0) = {y = 0, z = ±{x-xo), a;eK}. 

Existence of a one-dimensional manifold lTs/„(a:o, 0,0) that is tangent to Es/y_{xo,0,0) follows from The¬ 
orem 4.1 in [1]. We notice from (3.6) that y = 0 is an invariant reduction of the three-dimensional system. 
Therefore, we can set y = 0 and seek a parametrization of Wsf.^{xo, 0, 0) by working with z = '!/'(2:), where 
tpixo) = 0 and ^’'{xo) = ±1. Dynamics along Ws/u{xo, 0; 0) are given by the one-dimensional system 


m 


1 


X = -- 


xipfx). 


(3.16) 


From the three-dimensional system (3.6), we obtain a linear differential equation for nonzero ij}: 

df) m -\- 3 


dx m -I- 1 


'f’{x) ± : 


(3.17) 
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This equation completed with the initial condition tpixo) = 0 admits a unique solution 

1/ \ m+l 
V{x) = T—^— X 

Note that ip'{xo) = ±1, because lTs/„(a;o, 0,0) is tangent to 0, 0). Substituting 

1 p{x) = ±{x - Xq) + 0{{x - Xo)"^) 


1 - 

f-1 

2 - 

m + 1 


\XgJ 



(3.18) 


(3.19) 


into the differential equation (3.16) and truncating at the quadratic remainder term, we obtain the linear 


equation (3.15). Persistence of the linear dynamics on W);/„(xo, 0,0) with respect to the remainder term 
follows from analysis of the system (3.16) with the expansion (|3.19|). □ 


Proof of Theorem 1 1. ^[ For every Xq > 0, Proposition 3.1 states that there exists exactly one trajectory 
with y > 0 such that y(s) —>■ 0 as s —)■ -l-oo. This trajectory belongs to the center manifold Wc(a^o,0,0), 


whose dynamics satisfy the system (3.9). We recover the asymptotic behaviour (1.15) by eliminating y 


from the transformation (3.5) and setting x = Xg. 

For every xq > 0, Proposition 3.2 states that the point (a:o,0,0) is an intersection of the two- 
dimensional center manifold Wc(a;o,0,0) and the one-dimensional stable/unstable manifold Ws/u{xo,0,0) 
for the upper/lower sign. Therefore, the point can be reached in the direction s —>■ -|-oo along Wa{xg, 0, 0) 


but can not be reached along hF„(xo,0,0). This guarantees that the trajectory in Proposition 3.1 with 
y > 0 such that y(s) —)■ 0 as s —)■ -l-oo is unique for the lower sign. Therefore, for every Xg > 0, there 


exists a one-dimensional set of solutions of the differential equation (1.7) for the lower sign such that 
F7_(^) —y -)-oo as ^ —y -l-oo. 


On the other hand, the span of the trajectory in Proposition 3.1 and the trajectory in Proposition 3.2 
is a two-dimensional set that hosts all trajectories with y > 0 such that y{s) —)■ 0 as s —)■ -l-oo. Therefore, 


for every ccq > 0, there exists a two-dimensional set of solutions of the differential equation (1.7) for the 
upper sign such that i7+(C) —t +oo as ^ -boo. The rate of change along Wg{xg, 0, 0) is exponential in s 
and the rate of change along Wc(a::o,0,0) is algebraic in s. Therefore, solutions along the two-dimensional 
set still obey the asymptotic behaviour (1.15). □ 


4. Connection of invariant manifolds. 

on the existence of solutions H± 


Let us summarize the results of the previous two sections 
to the differential equations (1.7) on the semi-infinite line (7lj-,oo) that 


satisfy the properties (1.8)-(1.10). Such solutions are related to the trajectories of the dynamical systems 


(2.2) and (3.6), which depart from the equilibrium points where H± is zero and arrive to the equilibrium 


point where H± is infinite. In what follows, we will consider separately the two different systems for 
and H-. 


4.1. The system for H+ (t > 0). By Propositions 2.1 and |2.2[ for every nonzero A = A+, there is 
a unique trajectory of the dynamical system (2.2) in variables that departs from the equilibrium 

point (A+,0, 0) as T —>■ —oo and belongs to the domain m > 0. This trajectory is contained in the center 
manifold ITc(A+,0, 0) if > 0 and in the unstable manifold Wu(A+,0,0) if A+ < 0. 

By Proposition s |3.1| and |3.2[ for every xg > 0, there is a two-dimensional set of trajectories of the 
dynamical system ( |3.6| ) in variables {x,y,z) that reaches the equilibrium point (a:o,0,0) as s —)■ -l-oo 
and belongs to the domain y > 0. This trajectory is contained in the intersection between the center 
Wc(a;o,0,0) and stable lTs(a;o,0,0) manifolds. 

We shall now establish that the same trajectory departing from the equilibrium point (^+,0,0) in 


(2.2) arrives to the equilibrium point (xqiOjO) in (3.6). This trajectory determines a unique solution 
of the differential equation (1.7) with the upper sign satisfying properties (1.8)-(1.10). 


Lemma 4.1. Fix A^ € K\{0} and consider a one-parameter trajectory of the dynamical system (2.^ 
for the upper sign such that (^, u, w) —>■ (A_|., 0, 0) as t ^ —oo and u > 0. Then, there exists a tq € ffi (or 
Tg = -boo/ such that ^(r) — >■ -boo and u(t) — >■ -boo as t ^ Tg. 
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Proof. We introduce the energy-like quantity for the dynamical system (2.2) with the upper sign: 

1 


E{u,w) := -w^ - 

2 m - 


, m+1 


1 


m+2 


1 


TO -I- 2 


Computing the derivative of i? in r along a solution of system (2.2), we obtain 

dE 


1 


Cw . 


(4.1) 


(4.2) 


dr 2 

If A_|_ > 0, then ^(r) > A+ > 0, and i? is a strictly decreasing function of r as long as the solution 


to system (2.2) exists and w{t) ^ 0. By the representation (2.3) of VFc(^+,0, 0) in Proposition 2.1 we 
have w{t) > 0 for sufficiently large negative r. Now w{t) cannot vanish for any r, because if m = 0, then 
w = u"*(l -|- u) > 0, which contradicts positivity of w before vanishing. Therefore, E(w,u) decreases to 
—oo in finite or infinite time r. Because E{0, u) < E(w, u), we have u —>■ -l-oo if E{w, u) —>■ —oo. 

If A+ < 0, then E is an increasing function of r at least for sufficiently large negative r. By the 
representation (2.91 of W„(A+, 0,0) in Proposition |2.2[ we still have w{t) > 0 for sufficiently large negative 
T. Also, w{t) cannot vanish for any r by the same contradiction, since if re = 0, then w = u™(I -I- u) > 0. 
Therefore, ^(t) and m(t) are still increasing functions, and there is a finite tq € M such that ^(tq) = 0 
and f(ro) > 0. For r > tq, the energy method described above proves again that u{t) —>■ -l-oo in finite or 
infinite time t. 

We shall now prove that ^(t) —)■ -too as r —)■ tq. Since w > 0 for all r G (— oo, tq), the map t —>■ n is 
monotonically increasing, so that we can parameterize both ^ and w by u and consider the limit u —> -l-oo. 
From the last two equations of the system (2.2), we obtain 

dw , TO-I- I, 

w—=u (I-I-m)- ^w. 

du 2 

Since ^(r) > 0 for r close to tq and w(t) > 0 for all r G (—oo,to), we estimate 


(4.3) 


d 

du 


2 ^ 


< u™(l + u). 


Integrating this differential inequality, we obtain 

w^<C+ ^ 


,m+l 


,,m+2 


I 


TO -I- 2 


(4.4) 


where (7 > 0 is a constant of integration. Therefore, as rt —> oo, there exists a constant Wao > 0 such that 

m + 2 . . I I 

w < WooU^~ for sufficiently large u. From the first two equations of the system (2.2), we obtain 


d^ 

du 


u 

w 


> 


Wo 


Integrating this differential inequality, we obtain a lower bound for ^ given by 

2 


C>c- 


-U 2 


mwoo 


(4.5) 


(4.6) 


where C > 0 is another constant of integration. This lower bound yields ^ —>■ -l-oo as u ^ -l-oo. □ 


Corollary 4.2. The trajectory of Lemma f.l such that ^(r) — > -too and u{t) — > -|-oo corresponds 
to the trajectory of system [3.6] approaching the equilibrium state (a:o,0,0) for some Xq G [0, oo). Conse¬ 
quently, one can define a piecewise map 


D I — y Xq G 


(4.7) 
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Proof. We can use the transformation (3.51. Since u 
obtain 


z = wu 


< WooU 2 . 


00 , then y —)■ 0. From the bound (4.4), we 

(4.8) 


Therefore, z —> 0 as u —)■ oo. On the other hand, from the bound (4.6), we only obtain 


X = 2 > Cu ^, 

which is not sufficient to claim that x remains bounded as u —>■ oo. 


Using the transformation (3.5) and equations (4.3) and (4.5), we obtain the following system for x 
and z in the variable u: 


and 


dx 
du 

dz dx 1 
du du u^ 


1 

uz 


m - 


1 


1 

uz 


m + 3 


-uz 


(4.9) 


(4.10) 


It follows from the bound (4.8) and equation (4.10) that there is a positive constant C such that 


-—{z — x) > —Cu 2 , 
du 

for sufficiently large u. Since u~i is integrable asu —>■ -l-oo, then z —a; is bounded from below by a negative 
constant. Since z is bounded and approaches to zero as u —>■ oo, we finally obtain that x is bounded from 
above by a positive constant for all sufficiently large u. Finally, (xq, 0, 0) is an equilibrium state of system 
(3.6), therefore, the trajectory approaches (a;o,0,0) for some Xq € [0,oo). □ 


Remark 2. Integrating equation (f.9}, we obtain 


1 


u 2 du 


X = 


m + l 

U 2 


m + 1 

U 2 


UZ 


where c is an arbitrary constant. If we can prove that limu_ 

2 

lim xiu) = - 

u-foo aoo(m -I- 1) 


too uz{u) = Ooo > 0, then 
G (0,oo). 


This would indicate that the bound (4-.8) is not sharp. However, we only have numerical data supporting 
this claim for every G K. 


We will show numerically in §5 that both pieces of the map (4.7) are monotonically increasing for 
all G M and intersecting at xq = xq := y^2l{m -I- 1) for A+ = 0. Therefore, the piecewise map 

M\{0} 9 A^ i-A Xq G K"*" is in fact continuous at A_|_ = 0. The exact value xq corresponds to the exact 
solution (1.16) of the scalar differential equations (1.7). This exact solution is recovered with the following 


elementary result. 


Proposition 4.3. There exists an exact solution of the dynamical systems (2.2) and \ 3.6) given by 


(t^) ’ 

^ = ((Urf^) (t^) > 


G (-C 


(4.11) 
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and 


X{s) = XQ, 

yis) = % ^e(-6- 

z{s) = 


(4.12) 


respectively, where xq = whereas a and b are arbitrary positive parameters. 


Proof. Consider the system (3.6) and try the reduction x = xq, where xq is constant in s. Then, on 

cK 
2 


setting y = we obtain the following differential equations: 


, / 2 ^ 2 
y = -zy, z = y --—z . 


This system is compatible with the constraint y = ^^^^^Xgz if z' = —z^ and Xq = = Xq, so that 


y = XqZ. The general solution of z' = —z^ is z(s) = for a positive parameter b. The solution is 
defined for s > —b~^. 

Consider the system (2.2) and try the reduction = w = ii. Then, the system is compatible with the 
reduction if ^ = w"* and Therefore, 


e = 


m + 1\ 


which admits a general solution given by 


^(t) = 


2™(m + 1) 
(m - l)"*+i 


1 — ra 


n + l 
n-1 


where a is an arbitrary positive constant. The solution exists for t < a 
found from the above relations. □ 


Other components u and w are 


Remark 3. Note that in the exact solution (f. 7|), e(r)^0 = as r —)■ — oo. Therefore, the exact 
solution of Propositioncorresponds to the choice A_|_ = 0 in Lemma 4-1 


4.2. The system for H- {t < 0). By Propositions 2.1 and|2.2[ for every nonzero A = A-, there is 


a unique trajectory of the dynamical system (2.2) in variables {^,u,w) that departs from the equilibrium 


point (A_, 0, 0) as T —)■ —oo and belongs to the domain u > 0. This trajectory is contained in the unstable 
manifold W4(A_, 0, 0) if A_ > 0 and in the center manifold Wc{A_, 0, 0) if A_ < 0. 

By Propositio ns |3.1| and 3^ for every Xq > 0, there is a one-dimensional set of trajectories of the 
dynamical system ( |3.6[ ) in variables {x,y,z) that reaches the equilibrium point (a;o,0,0) as s —>■ -boo and 
belongs to the domain y > 0. This trajectory is contained in the center manifold Wc{xq, 0, 0). 

If we try an argument used in the proof of Lemma |4.1[ then it becomes clear that most of the 


trajectories departing from the equilibrium point [A- ,0,0) in system (2.2) will not arrive to the equilibrium 
point (xo,0, 0) in system (3.6) but instead reach the value u = 0 in a finite t S K. This indicates that, 
first, there are very few values of A-, for which the trajectories may reach infinite values for u, and 
second, the numerical method should not be based on the trajectories departing from the equilibrium 
point (A_, 0, 0) (such a shooting method was used previously in jS]). Instead, it may be better to look for 
the one-dimensional trajectory departing the equilibrium point (xo,0, 0) in system (3.6) in the negative 
direction of the time variable s. 


To illustrate the previous point, we show on figure 4.1 the trajectories of the system (2.2) with m = 3 


starting from the equilibrium point ( 2 l±, 0 , 0 ) along either center (for A+ > 0) or unstable (for A- > 0) 
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manifolds. The trajectories of the system for extend from small to infinite values of H+, see panel (a). 
Contrastingly, the trajectories of the system for H- turn back and return to small values of i?_, see panel 
(b). Note that the return time is significantly different between the first two and the last two trajectories. 
This indicates the presence of a particular value of A-, for which there exists a trajectory that extends 


from small to infinite values of H_, see (5.4) below. 




Fig. 4.1. Panels (a) and (b) show trajectories of the system with m = 3 for and H— respectively. In both 

cases, the trajectories start from the equilibrium point {A±, 0, 0) and depart along either the center or the unstable manifolds. 


In order to justify our numerical scheme, we shall prove that the trajectory departing the equilibrium 


w = 0 of system (2.2). 


point (aiojOiO) in system (3.6) in the negative s direction either intersects the plane u = 0 or the plane 


Lemma 4.4. Fix xq > 0 and consider a one-parameter trajectory of the dynamical system {S.tfjj for 
the lower sign such that {x,y,z) -A (xo,0, 0) as s ^ +oo and y > 0. Then, there exists an sq G M. (or 
Sq = —oo) such that 

(i) either 2 ;(so) = 0 and y{so) G (0, oo], 

(ii) or y{s) -A -boo as s ^ sq, whereas lim —G [0,oo). 

s—>-so y{s) 2 

Proof. For convenience, let us reverse the time variable, by transforming s —)■ — s, and rewrite system 


(3.6) with the lower sign in the negative direction of s: 


x' = '^xz - y, 
y' = zy, 

z' = y{l -y)- ^xz + ^z^. 


(4.13) 


By Proposition |3.1[ we have y > 0 and z > 0 for the trajectory departing from the equilibrium point 


(xo,0,0) (in negative s) along IFc(xo,0,0). From the second equation of system (4.13), y remains an 


increasing function of negative s as long as 2 ; remains positive. Therefore, we have an alternative: either 
z vanishes before y diverges or y diverges before z vanishes. 

The first choice of the alternative gives case (i). For the second choice, let us consider variables f and 
w given by (3.5) and parameterized by y in the limit y -G -boo (the map s i-G y is one-to-one and onto). 


By dividing the first and third equations in system (4.13) by the second equation, we obtain 


dx m -b 1 X 
dy 2 y 


dz 

dy 


1-y 


TO -b 1 X TO - 

2 y^~^ 


3 z 
y’ 
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By using variables ^ and w given by (3.5), we obtain 

_ 1 
dy y' 


m-\-2q 


dw 

dy 


1 — y TO + 1 ^ 


m+S^ 


y 


y"‘ 


(4.14) 


To show that the second choice of the alternative above gives case (ii), we will prove that w remains finite 
as y ^ + 00 . This is done by a contradiction. Assume that w —>■ +oo as y —?> +oo. Therefore, there exists 
Wo > 0 such that w > wq for all sufficiently large y. Then, from the first equation of system (4.14), we 
have for sufficiently large y, 


dy 


1 

“ y'^+'^wo' 


Because y 


— m—2 


is integrable at infinity, there is a finite positive ^oo such that |^| < ^oo for all sufficiently 


large y. Then, from the second equation of system (4.14), we obtain for sufficiently large y, 

y-l TO + 1 ^oo 


dw 

dy 


2/™+3wo 


yz 


Since both y~™‘~^ and y~^ are integrable at infinity, there is a finite positive Woo such that w < Wao for 
all sufficiently large y, contradicting the assumption that w —>■ +oo as y —?> +oo. Therefore, the case (ii) 
is proved. □ 


Corollary 4.5. The trajectory of system (3.6) departing from the equilibrium point (a:o,0, 0) with 
Xq > 0 in the negative direction of s intersects either the half-plane w = 0 and u > 0 in system 1(2.1 


in case (i) of Lemma 4-4 or the half-plane u = 0, w > 0 in system (2.2) in case (ii). Moreover, in the 
corresponding cases, 

(i) if lim u{s) > 0, then lim |^(s)| < oo 

S—^Sq S—^Sq 

(ii) if lim w{s) > 0, then lim |^(s)| < oo. 

S—^Sq S—^Sq 

Consequently, one can define two piecewise maps 


(i) 


9 xo i-A (^, m) S M X K+ and (ii) K'*- 5 y;) g x 


(4.15) 


Proof. In case (i), it is trivial to see that 2 ;(so) = 0 and y{so) € (0, + 00 ] corresponds to the half-plane 
w = 0 and u > 0. The first equation of system (|4.13 ) can be written for the variable f as follows: 



where the prime still denotes the derivative with respect to the time variable s in the negative direction 
of s. If y remains finite as s —)■ Sq (so that u(so) > 0), then ^(so) is bounded. 

In case (ii), it is also trivial to see that lims_>s y(s) = -|-oo and lim z(s)y(s) 2 ^ € [0, c») corresponds 

S— 


to the half-plane u = 0 and w > 0. If w remains nonzero in the limit y —> -l-oo, then, there exists Wg > 0 
such that w > wq for all sufficiently large y. Then, the same analysis as in Lemma |4.4| applies to the first 
equation of system (4.14) and shows that f remains finite as ?/ —>■ -l-oo. □ 


Unfortunately, we do not control the value of f at the intersection of the two piecewise maps (4.15). 
However, we will show numerically in ^that the piecewise maps (4.15) are typically connected at the 
points where u = w = 0 and ^ = A G M. In this case, a true solution iJ_ of the differential equation (1.7) 
with the lower sign satisfying properties (1.8)-(1.10) exist. 




























SELF-SIMILAR REVERSING INTERFACES 


17 


5. Numerical results. Let us describe a new numerical approach, based on the results of Lemmas 


4.1 and 4.4 that will be used to furnish meaningful solutions to the differential equations (1.7|, for 
and 77+. In §5.1| and §5.2| w e describe the numerical procedures for finding solutions for iJ_ and 77+ 
respectively. Finally, in §5.3| we summarize the results of the numerical experiments and compare them 
with the results found in Foster et. al. El- 


5.1. Solutions for 77_ (7 < 0). As discussed in 14.2 we wish to numerically construct a unique 
trajectory from infinite to hnite values of 77_. To do so, we integrate the system (3.6) from near the 


equilibrium point (a;o,0, 0) in the far-held towards the equilibrium point (A_,0,0) of the system (2.2) in 
the near-held. The numerical procedure is carried out as follows: 

1. Select a value of xq > 0. Ideally, one would like to begin by using this choice of xq to specify 


unique initial values for {x,y,z) = (a;o:0,0), and then numerically integrating the system (3.6) backward 
in the ‘time’ variable s. Equivalently, one could integrate the system ( 4.13| ) forwards in time. However, 
since {xq, 0,0) is an equilibrium point, it is not possible to escape (xq, 0,0) in a hnite time. Thus, in order 
to ensure that any numerical integration scheme can depart from near the equilibrium point along the 
center manifold, Wc(a:o,0,0), it is necessary to take a ‘small step’, say 5^1, away from (a:o,0, 0) using 


the relevant asymptotic behaviour. Using (3.7) and (3.8) we hnd that a trajectory on the center manifold 
114(a:o,0,0) has the local behaviour 


X = Xo+ - ( 

y=^xoS + 0(S^), 
z = S, 


m+1 

2 


fx^,)s + 0(S% 


(5.1) 


for small positive values of S. Having selected values for both xg and S, the behaviours (5.1) may be used 
to specify unique (pseudo-)initial values for (x, y, z) and to begin the numerical integration of the system 


(3.6) in the direction of decreasing s. In this study, numerical integration of the system (3.6) was carried 


out using the ode45 routine in MATLAB with the default settings, except AbsTol and RelTol which were 
both set to have a value of 10“^°. Selecting an appropriate value of 5 is a somewhat ad-hoc procedure: 
there is trade-off between taking <5 too small, which renders it difficult for the numerical integration to 
escape the neighbourhood of the equilibrium point (leading to poor accuracy of the integration), and taking 


5 too large which could result in low accuracy of the asymptotic expansion (5.1). However, we found that 
choosing (5 € (10“^, 10“^) gave good results over the ranges of parameters we studied. Robustness of the 
results with respect to changes in: (i) the choice of 5, and; (ii) the number of terms in the asymptotic 
expansion were verified. 

2. The result of Corollary |4.5| asserts that all such trajectories will ultimately — in either a finite or 


infinite time — intersect with either the plane w = 0 or u = 0, see the maps (4.15). To ensure accurate 
numerical integration of the system in the near-held, in variables it is necessary to ‘switch’ 


from integrating the far-held system (3.6) to the near-held system \2.2\ backward in time. The choice of 
conditions under which this switch should occur is, again, somewhat arbitrary. However, as long as the 
values of {x, y, z) — and hence the values of (^, u, w) — are all hnit e an d no n-ze ro, this switching is valid 
at any point. In this study, we chose to switch from integrating (3.6) to (2.2) when = 20 

(or equivalently when ^ = 20). However, we verihed that our results were robust to changes in the 
choice of switching conditions. This switching procedure can be readily automated within MATLAB using 
the Events function to autonomously: (i) stop the integration of the system (3.6) when specihed the 
conditions are satished; (ii) read-off the hnal values of {x,y,z); (hi) transform these to corresponding 
values for {^,u,w) using the change of variables (3.5), and; (iv) begin the integration of (2.2) backwards 
in time from the appropriate initial data. 


3. The integration of the system (2.2) is then continued backward in the time variable t until either 
w = 0 or u = 0. Again, we used the Events function to autonomously detect when either of these events 
occurred and to stop the integration. It is noteworthy that we found it helpful to use odelBs to integrate 
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the near-held system — again, the default settings were used with the exception of AbsTol and RelTol 
which we both set to 10“^°. Although odel5s is typically slower than ode45, it is considerably more 
appropriate to deal with integrating systems of equations that exhibit apparent stiffness. This apparent 
stiffness, manifested as rapid changes in the direction of the trajectory in variables arises as 

an artifact of the inhnite time required to reach the equilibrium point (A_,0,0). Thus, if a trajectory 
approaches very close to the equilibrium point it appears to be rapidly rejected from that neighbourhood. 
When the integration is terminated, we record the following pieces of data: (i) the selected value of xq; 
(ii) whether the trajectory reached m = 0 or w = 0, and; (hi) the value of either (^, u) or ( ^, w) at the 
termination point. These data define a point on one of the two piecewise maps defined in (4.15). It is 


by computing a large number of trajectories, each emanating from different values of xq, that we are able 
to trace out the forms of these maps. 

The non-trivial solutions for iA_ that we seek correspond to trajectories emanating from particular 
equilibrium points in the far-held, say (a:g,0,0), that reach the near-held equilibrium point, (A_,0,0) for 
some hnite A_ ^ 0. In addition to th ese n on-triv ial solutions, we recall that for all values of m > I there 


exists a trivial solution given by (4.11) and (4.12) that emanates from Xq = xq = and reaches 

A_ — 0, as discussed in Proposition 4.3 Some representative results for m = 2,3,4 and 5 are shown in 


hgure [5T| In these plots, a suitable non-trivial solution for iJ_ is found by identifying a value of xq = Xq 
for which the value of (^,u) = (A_,0) — or (^,w) = (A_,0) — at the termination point. 

In hgure we show some representative computations highlighting the differences between the near- 
held behaviour of trajectories local to a true solution with A_ < 0 and A_ > 0. As is e viden ced by hgure 
5.1 close to a trajectory with A_ > 0 the piecewise continuous maps dehned in (4.15) are smooth, 


whereas close to a trajectory with A_ < 0 the maps exhibit rapid changes and (vertical) cusp-like features. 
As a result, determining negative value(s) of A_ is more challenging — despite resolving xq to machine 
precision (approximately 10“^^ in the standard IEEE double precision), it is not possible to approximate 
the value of Xq sufhciently well that an accurate estimate of A_(< 0) can be determined. In these cases, 
we therefore found it necessary to implement one additional stage in the numerical scheme as follows. 

Once Xq had been determined up to machine precision, and two ‘limiting’ trajectories had been 
identihed (one emanating from a;g ± (5 and terminating at u = 0, and the other emanating from ccg <5 and 
terminating at w = 0, where 5 <C 1) the expected near-held linear asymptotic behaviour of u(^), according 
to the expansion (1.13) for the solution H- in Theorem is clearly visible. This linear behaviour can 
then be extrapolated, in the direction of decreasing until it intersects the ^-axis. This intersection point 
is, to a good approximation, the value of A_(< 0) corresponding to the trajectory emanating from Xq. 
Panel (a) of hgure 5.2 gives an example of the linear extrapolation procedure described above. 

Using the procedure described above, for all values of m > 1, we recover the trivial solution discussed 
in Proposition |4.3[ In addition, we see that in the case m = 3 there is only one suitable non-trivial solution 
with the following data: 


m = 3 : 

* ^ 
Xq ; 

« 0.767, 

A_ 

« 0.129. 

For m = 4 and m = 5, similar results are observed with 
follows: 

one 

trivial and only one non-trivial solution 

m = 4 : 

* ^ 
Xq ; 

« 1.165, 

A_ 

« 0.386, 

and 





m = 5 : 

* ^ 
Xq ; 

« 1.666, 

A_ 

« 0.501. 

Contrastingly, in the case m = 2 we hnd that three suitable non-trivial solutions exist with the data: 

1 

r 

:! 0.338, 

A_ 

« -2.804, 

m = 2 : 1 

X^ R 

:! 0.137, 

A_ 

« -0.932, 

1 

1 ^0 « 

:! 0.0592, 

A. 

. « -0.546. 
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(a) (b) 




Xo 


Xo 


(c) 



Xo 


(d) 



Xo 


Fig. 5.1. Panels (a)-(d) show plots of the piecewise maps defined in l[4-b5]) for m = 2,3, 4 and 5 respectively. In 
all cases the blue, red and black curves show the value of w at u = 0, the value of u at w = 0 and the value of ^ at the 
termination point respectively. The dashed vertical line indicates the value of xq = xq corresponding to the exact solution 
Hi.16^. The crosses on panel (a) mark the data points extracted using the extrapolation procedure discussed in the text. 


Notably all values of A- for m = 2 are negative, whereas for m = 3,4 and 5 they are positive. 

In addition to the detailed results for m = 2, 3, 4 and 5 shown in figure [O] we also show in figure 
5.3 the values of Xq determined for all values of m from m = 1 to m = 8. For the primary red and blue 
branches, emanating from m = 3 along the black branch, we show the corresponding values of A- in figure 
|5.4[ Intriguingly, the numerical results indicate that in addition to the exact solution - which is valid for 
all TO > 1 — there are a whole host of additional solutions, some with A_ > 0, and others with A_ < 0. 
In particular, there is at least one additional trajectory corresponding to a suitable solution for H- with 
A- > 0 for all values of to > 2.978. Further, for all to < 3 there exists at least one additional solution for 
i7_, although, in this case for a value of A- < 0. Another noteworthy feature of the plots shown in panels 
(a)-(d) of figure 5.3 is that at each value of to = {2N — I) for TV S N additional branches of solutions 


depart from the branch along which xq = xq and A_ = 0. The underlying reason for this structure is as 
yet not understood. 

5.2. Solutions for H+ {t > 0). Having successfully found suitable solutions for iJ_, we now proceed 
to compute suitable solutions for iJ+. As discussed in §4.I[ we can numerically construct a unique tra¬ 


jectory from small to infinite values of 77+. To do so, we integrate the system (2.2) from the equilibrium 
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Fig. 5.2. Panel (a) shows some representative trajectories emanating from Xq ^ 0.338 for m = 2. More precisely 
the red and blue trajectories begin at Xq ± S for 5 = 10“^, 10“^, 10“®, 10“®, 10“^^ and 10“^^. The black line shows the 
artificially extrapolated linear behaviour. Panel (b) shows some representative trajectories emanating from Xq ^ 1.165 for 
m = 4. In this case, the trajectories begin at Xq ± S for 5 = 10“^, 10“^ and 10~®. Notably, in the latter case, despite only 
resolving Xq to 3 significant digits, a good estimate of A- has already been obtained. 


point (^+, 0,0) in the near-held towards the equilibrium point (xq, 0, 0) of the system (3.6 1 in the far-held. 
The numerical procedure is carried out as follows: 

1. Select a value of A+ e K\{0}. Since (A+, 0, 0) is an equilibrium point, it is not possible to escape 
(A+jO, 0) in hnite time. We therefore begin integration of the system (2.2) by taking a small step, say 
e, away from (yl+,0, 0) using the relevant asymptotic behaviour. Using (2.3) and (2.4) for A+ > 0, we 
hnd that a trajectory exiting the equilibrium point along the center manifold, Wc(^+,0, 0), has the local 
asymptotic behaviour 



for A+ > 0, 


+ 0 { 






(5.2) 


for a small positive value of e. By contrast, using (2.9) and ( 2.10[ ) for < 0, we hnd that a trajectory 
along the unstable manifold Wu{A+,0, 0) has the local asymptotic behaviour 


u= (m^\A+\m)^ e 


0{e 


for A+ < 0, 


(5.3) 


0{e) 


for a small positive value of e. Having selected values for both A^ and e, either (5.2) or (5.3) dehne 


unique (pseudo-)initial conditions to begin integrating the system (2.2) in the direction of increasing time 
T towards the far-held. 

2. We proved in Corollary |4. 2 [ that the ultimate fate of all such trajectories, in variables {x,y,z), 
is approachi ng t he equilibrium state (aio,0, 0) for some Xq G [0, oo). Thus, by continuing integration of 
the system (2.2) to some large value of r, denoted by say Too, and reading off the value « xq 

at T = Too, we can obtain an arbitrarily accurate approximation of the corresponding value of xq that is 
obtained in the far-held — a higher degree of accuracy can be achieved by simply increasing the value of 
Too- For this purpose we found ode45 with the majority of the default setting to be sufficiently robust. To 
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(b) 


(c) 


id) 





Fig. 5.3. The variation of Xq versus m. The red, blue and black curves indicates values of Xq that define trajectories 
terminating at the near-field equilibrium point with A— > 0, A- < 0, and A— = 0 respectively. Panels (b)-(d) show 
zoomed-in regions from panel (a) near m = 3, m = 5, and m = 7 respectively. 


ensure high numerical accuracy, at the cost of a relatively small increase in computation time, both AbsTol 
and RelTol were decreased to 10“^°. In contrast to the case for solutions i7_, we found it unnecessary to 
‘switch’ from integrating the near-held system (2.21 to the far-held system (3.61. Typically, we found that 
taking Too & (10^, 10®) gave an approximation of xq correct to 8 signihcant digits. 

Carrying out this procedure for a variety of choice s of A+ we ar e ab le to trace out the form of the 
piecewise map between and xq dehned earlier in (4.71. In hgurewe show this map for m = 2, 3 
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Fig. 5.4. Panel (a) shows the variation of A— along the red and blue curves emanating from the black curve near 
m = 3. Panel (b) shows the same plot zoomed in on positive values of A—. 




Fig. 5.5. Panel (a): Plots of the variation of xq with A+ for various different values of m = 2, 3 and 4. Panel (b): 
Plots of the trajectories emanating from A+ = 0.1, 0.2, 0.3, 0.4 and 0.5 for m = 3. The constant to which these trajectories 
tend in the far-field is selected to be the corresponding value of xq- 


and 4 (see panel (a)), as well as some representative trajectories of the system (2.2) for m = 3 (see panel 


(b)). In addition to the results shown, other computations for different values of m were also carried out 
and it appears generic that Xg is a monotonically increasing function of A^. Crucially, it appears that 
range of the map (4.7) is the entire semi-axis K+ for xg- 


5.3. Summary of numerical results. We have demonstrated that: (i) for each value of m > 1 
there exists at least one value of xg = Xg (different from the trivial case Xg = xq) that defines a trajectory 
emanating from (xg, 0, 0) and terminating at a (A_, 0,0), and thus a suitable solution for 7J_, and; (ii) for 
every value of S K there exists a unique corresponding value of Xg, thereby defining an infinite family 
of suitable solutions for H^. The one remaining step is therefore to invoke the matching condition (1.11). 


This condition is equivalent to requiring that the far-held behaviour of is characterized by xg = Xg. 
Thus, given a solution for i7_, the matching condition ( 1.11 ) specihes a unique choice of a;o = Xg, a unique 
A+, and thus a unique solution for i7+, thereby closing the problem. 
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The solutions found here for m = 3 and 4 show a qualitative, although not quantitative, agreement 
with those reported in [B]. Here, we found that 


m = 3: 0.129, H+« 0.154 and Xq « 0.767 (5.4) 

and 

TO = 4 : « 0.386, A+ « 0.794 and « 1.165. (5.5) 

In Foster et. al. [B], they claimed that for to = 3: A_ « 0.144, H_|_ « 0.0958, and Xq « 0.765, whereas for 
TO = 4: A- « 0.386, A^ « 0.341, and Xg « 0.980. Additionally, they claim another suitable solution for 
TO = 2: A_ « 0.00135, A_|_ « 0.0102, and Xq « 0.817. In contrast, here we found that no such solution 
with A_ > 0 exists. Notably this value of a;o reported in for m = 2 is very close to the value of 
Xq = y/2j{m + 1). For to = 2, using our numerical approach we have been able to identify three other 
solutions with A_ < 0, namely: 


1 

r 

» -2.804, 

A+r 

i -4.322, 

Xq « 

0.338, 


TO = 2 : 1 


» -0.932, 

A+s 

i -30.625, 

Xq - 

a 0.137, 

(5.6) 

1 


a -0.546, 

A+« 

:! -166.623, 

Xq 

« 0.0592. 



We believe that the origin of these discrepancies is due to the low accuracy of the numerical scheme used in 
[B]. Indeed, in [B], the solutions for i7_ were computed by identifying the value of A_ which characterizes 
solutions in the near-held that extend into far-held with the requisite behaviour, as in panel (a) of hgure 
|4.1[ Solutions for were computed by hnding the value of A^ inferred (via shooting from the far-held 


toward the near-held) by invoking the matching condition in the far-held as in panel (b) on hgure 5.5 
Both numerical methods used in [B] are ill-posed. Here, we pose the numerical problem as a shooting 
scheme for uniquely dehned piecewise scalar functions, i.e. the maps dehned in (4.71 and (4.15). We 
therefore believe that the results obtained here are more reliable than those in [BI. 




Fig. 6.1. Representative plots of h(x,t) at 10 equally spaced values of t between —1 and 1, and for t = ±10”^ 
to demonstrate the continuity of h[x,t) across t = 0. Solid and dashed curves show the solution for t < 0 and t > 0 
respectively. Panel (a) shows the anti-reversing dynamics for m = 2, A— —2.804, A+ ps —4.322 and Xq ^ 0.338. Panel 
(b) shows the reversing dynamics for m = 4, A— ps 0.386, A+ 5^ 0.794 and Xq ps 1.165. 


6. Conclusion. This work has focused on constructing local (in both space and time) self-similar 
reversing and anti-reversing solutions to the nonlinear diffusion equation with TO > 1 and n = 0. We 
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have demonstrated how the dynamical theory combined with the numerical scheme can be used to furnish 
suitable solutions to the differential equations (1.7) for i7_ and H+. Via the self-similar reductions (1.6), 


the solutions to these differential equations can be transformed into physically meaningful solutions for 


h{x,t) to the nonlinear diffusion equation (1.1) with m > 1 and n = 0, which is completed with the no 


flux boundary condition (1.3) and the condition the = 0. In this final section we shall discuss the 

connection between the self-similar solutions found here and the dynamics of the model (1.1). 


It is well-known, and can be readily verified, that the nonlinear diffusion equation with absorption 
(1.1) admits two different travelling wave solutions. Seeking such solutions near a left interface x = £{t) 
that have the form h{x,t) ~ A{t){x — for some function A{t) and constant a, we find that 


as x\£{t), for £ < 0, 

h ^ {£)~^{x — £{t)) as x\£{t), for £ > 0. 


( 6 . 1 ) 

( 6 . 2 ) 


The former, is an advancing wave local to a left interface whose motion is driven by diffusion, whereas 
the latter is a receding wave driven by absorption. This study has therefore elucidated the process by 


which the wave (6.1) becomes (6.2), giving rise to a reversing interface, or vice versa, giving rise to an 


anti-reversing interface. 


In (5.1 we used the result of Lemma 4.4 to numerically construct suitable solutions for i7_. For each 
value of TO > 1 we identified at least one suitable solution, defined by a pair of values of A- and Xq, 
in addition to the exact solution dlT^ — in the original time and space variables, this exact solution 
corresponds to a steady solution for h{x, t) and thus does not constitute a reversing nor an anti-reversing 
solution. In j ]5.2[ we used Lemma |4.1| to formulate a numerical scheme for constructing solutions for 
defined by pairs of values of A+ and Xq. We showed that the map (5.5) is one-to-one and its range is 


the entire semi-axis K+ for xq. Importantly, for each value of to > 1, we found that: (i) if A^ < 0 then 
Xo < XQ, and (ii) if A^ > 0 then Xq > xq where xq = -I- 1). The final stage in constructing 

solutions for h(x,t) is to invoke the matching condition (1.11) that ensures continuity of h{x,t) across 
t = 0. Owing to the aforementioned properties of (5.5) we are forced to reject any solution for H+ that is 


defined by a trajectory with A- < 0 and Xq > xq or A- > 0 and Xq < xq on the basis that it necessarily 
cannot match to solution for iJ_. In summary we have found that for 1 < to < 3 up to 5 different 
solutions are available with A-,A^ < 0. For to > 2.97 there is at least one solution with A-, A^ > 0. For 
7 < TO < 7.75 an additional branch of solutions with A_,A^ < 0 emerges, whereas for to > 6.42, there is 
another branch of solutions with A_, A+ > 0. For to > 8 it seems quite possible that yet more branches 
of solutions will emerge. 

There is a distinct difference between the interpretation of solutions with A- , A+ > 0 in terms of the 


original model (1.1) compared to those with A-,A^ < 0. The former, correspond to a reversing solution 
where the left interface advances for t < 0, with the behaviour (6.1), and then subsequently recedes for 


t > 0 with the behaviour (6.2). Contrastingly, the latter corresponds to an anti-reversing solution where 


the interface recedes for t < 0, with the form (6.2), and then advances according to the form (6.1) for 


t > 0. One representative local solution for h{x,t) for both types of behaviour — one reversing and one 
anti-reversing — are shown in figure |6.1[ 

Some natural open questions raised by this study are: (i) whether any self-similar solutions with 
non-monotone profiles (in ^) exist — i.e. solutions that do not satisfy (1.9); (ii) whether the self-similar 
solutions identified here are stable in the context of the model (O’ and; (iii) if more than one reversing 
or anti-reversing solution is stable for a particular value of to, what is the mechanism for selecting the 
appropriate self-similar solution at a particular reversing or anti-reversing event. 
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